Topological defects in the crystalline state of one-component plasmas of non-uniform 
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We study the ground state properties of classical Coulomb charges interacting with a 1/r potential 
moving on a plane but confined either by a circular hard wall boundary or by a harmonic potential. 
The charge density in the continuum limit is determined analytically and is non-uniform. Because of 
the non-uniform density there are both disclinations and dislocations present and their distribution 
across the system is calculated and shown to be in agreement with numerical studies of the ground 
state (or at least low-energy states) of A'^ charges, where values of A'^ up to 5000 have been studied. 
A consequence of these defects is that although the charges locally form into a triangular lattice 
structure, the lattice lines acquire a marked curvature. A study is made of conformal crystals to 
illuminate the origin of this curvature. The scaling of various terms which contribute to the overall 
energy of the system of charges viz, the continuum electrostatic energy, correlation energy, surface 
energy (and so on) as a function of the number of particles A'^ is determined. "Magic number" 
clusters are those at special values of A^ whose energies take them below the energy estimated from 
the scaling forms and are identified with charge arrangements of high symmetry. 

PACS numbers: 61.72.Bb, 61.72.Mm, 61.72.Lk 
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I. INTRODUCTION 

Classical charges moving on a plane and repelling each 
other via a Coulomb 1/r potential have a ground state 
which is a triangular lattice (for appropriately chosen pe- 
riodic boundary conditions). This is a two-dimensional 
(2D) example of Wigner crystallization and is known 
to occur in diverse areas of physics, such as electrons 
trapped in surface states of liquid helium colloidal 
suspensions and quantum dots. In this paper we exam- 
ine situations when the charges do not have a uniform 
density across the system. This occurs for example when 
the charges are confined by a circular hard wall or when 
they are confined by a harmonic potential. 

From a geometrical view point a perfect crystal lattice 
is a periodically repeating arrangement of identical struc- 
tural cells, which fit together without gaps or overlap. 
The question we seek to answer, for a classical Wigner 
crystal of non-uniform density, is: how much of what we 
mean by "crystal lattice" still applies to the resulting 
structure? On the one hand, one expects the structure 
of the lattice to be locally triangular; so that each lattice 
site has six nearest neighbors, since this is the optimal 
energy arrangement for uniform density. On the other 
hand, due to the changing density not all the symmetries 
of the triangular space group, such as the translational 
and rotational invariances can continue to apply. We 
seek to understand how this conflict is resolved. We are 
particularly interested in knowing if the resulting struc- 
ture can be understood within the framework of elasticity 
theory, and if so, what role is played by plastic deforma- 
tions, such as dislocations and disclinations. In addition, 
using a continuum model, we shall attempt to quantify 
the scaling with N of various phenomena associated with 
the structure of the cluster of charges, such as the energy, 
correlation energy (see below for its definition), surface 



energy and so on. This will in turn allow us to develop 
a link between symmetry and the energy of the lattice; 
we expect states with a high degree of symmetry to have 
a particularly low energy. Such states are known in the 
literature as "magic number" states. 

An experimental realization of a system of 2D charges 
confined by a hard wall might be electrons trapped in 
surface states of liquid ^He [T] . The hard wall potential 
could be effected by a circular boundary made from an 
electrical insulator. A harmonic confining potential also 
produces a non-uniform density across the system and 
has been studied by Koulakov and Shklovskii, see 0] and 
Their study is relevant to the properties of quantum 
dots. It was found that the density of the charges was 
greatest at the center and diminished, upon approach- 
ing the edge (whereas with hard- wall confinement, the 
density rises from the center of the disc towards the wall 
because of the repulsion between the charges). Their sim- 
ulations showed that although the charge density was not 
uniform, nevertheless the lattice was locally triangular. 
They showed that this is possible because the changing 
lattice density is accompanied by plastic deformations. 

It is a consequence of topology that both the hard- 
wall and the harmonic potential systems must contain an 
excess of 6 positive disclinations or pentagonal regions. 
Furthermore, the changing density introduces disclina- 
tions throughout the cluster. Disclinations also occur at 
the edges of the clusters to allow the lattice to adapt to 
the imposed circular structure. In addition to disclina- 
tions, the cluster includes dislocations (a tightly bound 
five-seven coordinated disclination pair). Unlike discli- 
nations the total number of dislocations is not fixed by 
topology; dislocations are present to reduce the strain en- 
ergy in the crystal which is induced by the circular edge 
and the disclinations. Dislocations are present in large 
numbers near the lattice edge, where they form a cloud 



2 



around any disclinations and help reduce the large elastic 
stress induced by the latter. However, it was found by 
Koulakov and Shklovskii that disclinations and disloca- 
tions are to be found also in the lattice interior and act 
as a mechanism for reducing the strain energy there. 

A number of other authors have carried out numerical 
simulations in situations where the classical 2D Wigner 
crystal has a non-uniform density. Bedanov and Peeters 
[3| considered particles interacting via the pure Coulomb 
potential and confined either by a hard wall or a parabolic 
potential. Simulations on the hard wall problem have 
also been carried out by Kong et al Q; in addition to 
the Coulomb interaction they also consider the dipole 
and Yukawa potentials. Ying-Ju Lin and Lin I 0] have 
studied a number of systems with different interaction 
and confining potentials. Our work confirms and ex- 
tends these earlier numerical studies, but also includes 
an account of our attempts to understand the numerical 
results. 

The paper is organized as follows. In Sec. II, we 
present the systems we study and our numerical ap- 
proach. In Sec. Ill we expand on the work by Koulakov 
and Shklovskii and derive the continuum approximation 
for the cluster of charges. We can then give an analyt- 
ical calculation for the density of charges in the contin- 
uum limit, and for the density of dislocations and discli- 
nations. We give a series expansion which we believe 
describes the various contributions to the ground state 
energy of the cluster. These results are compared to nu- 
merical experiments in Sec. IV. The discussion is in Sec. 
V, especially of the striking lattice curvature effect visi- 
ble in our studies of large clusters which is compared to 
that seen in conformal crystals. 

II. NUMERICAL APPROACH 

The energy of a cluster of N charges confined to a disk 
of radius R, by a hard-wall potential is given by 

N N 

where 

, X / for < i? 
~ \ oo for > R. 

The energy of a cluster of N charges in a parabolic con- 
fining potential is given by 

N N 

where we set A=l/2. 

Finding the global minimum for a function such as 
or is a very difficult task. The number of metastable 
states proliferate exponentially with N; consequentially 



the global minimum is obscured by a vast number of lo- 
cal minima with energies close to that of the global mini- 
mum. There exist a number of heuristic methods for such 
problems. Although there is no guarantee of finding the 
global minimum, it is possible to find states close to it. 

We found that for the hard-wall system the standard 
Metropolis simulated annealing algorithm to be more ef- 
fective than a conjugate gradient algorithm [7]. For a sys- 
tem with N charges the simulated annealing algorithm 
was run with typically x (5 x 10^) Monte Carlo steps. 
The temperature of the simulation was decreased linearly. 
The average displacement of the charges at each temper- 
ature step was chosen by an automatic process to give 
an acceptance probability of 0.5 ± 0.01. Promising states 
were reheated and annealed repeatedly to iron out as 
many defects as possible. Finally the results were put 
through a conjugate gradient algorithm to remove any 
residual strains. 

For the harmonic (parabolic) potential case, we found 
the conjugate gradient algorithm to be as effective as sim- 
ulated annealing. We used the former method for this 
system as it ran faster. Results were generated starting 
from an initial random configuration, which had a ra- 
dial density profile matching the continuum limit density 
given below in Eq. p^ . 

III. THE CONTINUUM LIMIT 

In the following we develop the continuum model of 
the two systems. For the case of parabolic confinement, 
many of the important results have already been derived 
by Koulakov and Shklovskii, see and so where ap- 
propriate we shall simply quote the relevant result. The 
bulk of the material in this section is concerned with the 
system with a hard-wall confining potential. 

In the following, for the system with the hard-wall con- 
fining potential, we use a variational approach to derive 
the (non- uniform) charge density in the continuum limit. 
Next we demonstrate that as a consequence of the non- 
uniform density the system interior will contain topolog- 
ical defects, where the density of these defects depends 
on the rate of change of the density. Finally, for the 
system with a hard-wall confining potential, a series is 
developed which includes the contributions to the energy 
of the cluster which scale smoothly with system size. 

A. Charge Density 

For charges confined by a hard-wall potential, the en- 
ergy expression in Eq. ([1]) can be approximated by the 
integrals over the disc r < R, 

The continuum approximation treats the density p^(r) 
as a smooth function rather than the sum of delta func- 



3 



tions 



N 



(4) 



where is the position of the i*'' charge. One then 
minimizes the energy of the cluster, with respect to the 
smooth function p^(r), subject to the constraint that the 
number of particles 



N : 



(5) 



is constant. Introducing the Lagrange multiplier fi the 
constrained equation is 



E 



1 /• ,2 P"ir) 

2 I '"ir-r'l 



A variation in the energy is given by 

6E = E[p"{r) + Sp"iv)]-E[p"{r)], (6) 

where 5p^{r) represents a small change in the charge 
density. Keeping only terms up to first order, Eq. ([6]) 
gives 



6E 



5p"{v') 



p"{v) 



(7) 



To make the functional derivative stationary we require 
'd\J^-p = 0. (8) 



To solve this integral equation it is convenient to write 
the integral in Eq. ([8|) in terms of radial and angular 
variables: 



d9 



p^ {r)rdr I ^ 

^0 (r^ + r' — 2rr'cos( 

R irp^ir) ( 2^/V?\ 

dr —K , 

r + r' \ r + r' / 



where K{k) is an elliptical integral of the first kind. This 
is a Fredholm integral equation of the first kind which 
has the radially symmetric solution Q 



P"{r) 



1 



(9) 



To determine the Lagrange multiplier we substitute Eq. 
((9]) into Eq. ([5]), which gives p = Ntt/2R, and we finally 
have 



N 



1 



27ri?2 



(10) 
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Incidentally, this result is the density profile obtained if 
a hemispherical shell of charge is projected onto a plane. 



It is important to note the physical meaning of the 
Lagrange multiplier p; it is the electric potential at any 
point inside the disk when the charge is distributed ac- 
cording to Eq. (|10p . This can be seen by referring back 
to the original formulation of the problem as given by 
Eq. dH). 

The number of charges within a distance r from the 
center of the disc is given by integrating Eq. (jlOp over 
the region r' < r: 

N"{r)^N(^l-^l-[^yy (11) 

The charge density in the continuum limit for a cluster 
of charges in a parabolic confining potential is [3| 

p^{r)^Po\ll^{^)\ (12) 

where 

The total number of charges within a distance r from the 
center of the disc is 



(13) 



N^ir) = iV I 1 - ( 1 - ^- 



B. Density of Defects 

Volterra dislocations are plastic imperfections charac- 
teristic of a deformed solid . For a 2D lattice the only 
relevant Volterra dislocations are the edge dislocations 
and the wedge disclination; these we discuss in turn be- 
low. 

A dislocation in a perfect lattice can be created by 
making a cut in the lattice, translating the cut edges with 
respect to each other and inserting/removing material. 
This process can also be viewed as the insertion/removal 
of a half plane of atoms, it is characterized by a discrete 
Burgers vector B which measures the amount by which 
the Burgers circuit around the dislocation fails to close 

A disclination can be created in a perfect lattice by 
making a cut, rotating the cut edges with respect to each 
other (this then defines a disclination axis) and insert- 
ing/removing a wedge of material. After the wedge is 
inserted/removed the whole construction is welded to- 
gether and allowed to relax. The closure failure of a 
Burgers circuit around a disclination is given by 



B = O X r. 



(14) 



where r is the distance from the disclination axis Jl, and 
|r2| is the wedge angle [ll|- In a real lattice the wedge 
angle is quantized by the lattice symmetry, thus in a 
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triangular lattice matter is inserted into the lattice if 
|r2| = +27r/6 and removed if |J7| = — 27r/6, which corre- 
sponds to a heptagon or a pentagon in the crystal lattice 
respectively. 

It is important to know the relationship between dislo- 
cations and disclinations. A dislocation is a tightly bound 
pair of disclinations of the opposite sign. A disclination 
can be decomposed into a series of dislocations, each of 
which have the same sign. For an illustration of both 



these points see [11 



Disclinations can be present in the ground state of a 
system for two reasons; either because they are demanded 
by topology (a consequence of Euler's theorem), or be- 
cause the lattice density is non-uniform. 

The first point has been covered in depth elsewhere 
Q, it will suffice to say that a Delaunay triangulation 
of a lattice will produce a unique planar graph. Euler's 
theorem states that for any such graph in flat space the 
following relationship holds between the number of ver- 
tices V, edges e and faces f 



v + f 



1. 



As a consequence of applying Euler's theorem to a trian- 
gular lattice we can assign a topological charge to each 
lattice site (or vertex), the sign and magnitude of the 
charge depends on by how much the coordination number 
(i.e. the number of nearest neighbors the site has) differs 
from 6. For example a pentagon has a topological charge 
+1 while a square has a charge of +2. Similarly a hep- 
tagon has a topological charge of -1 while an octagon has 
-2. Obviously a hexagon is topologically neutral. The to- 
tal topological charge for any cluster is always conserved 
and must always be equal to -|-6. 

We now show that a change in the density of the lattice 
leads to a closure failure of the Burgers circuit. This will 
allow us to calculate the Burgers vector density, from 
which the density of dislocations and disclinations can 
be determined. 

Consider a lattice with smoothly varying density and 
let us suppose that the increase in density depends only 
on the radial distance r. Drawing a square of dimensions 
Ar around such a region of lattice, see Fig. [U we see that 
the number of lattice rows crossing the side pq is given 
by 

L{r) 



air)' 



and the number crossing the side sr is given by 



L{r + Ar) 



Ar 



a{r + Ar) ' 



where a(r), the lattice spacing, is a function of density 
and in a triangular lattice is given by 



a(r) 



p{r)V^' 



(15) 



Ar 



a(r) 



|B(r)| 



[a(r+ Ar) 



Ar 



FIG. 1: Consider a region of crystal in which the lattice den- 
sity is increasing in the r direction. It can be seen that by 
drawing a square Burgers circuit (of dimension Ar) around 
the region that there are more lattice lines crossing the side 
sr than crossing the side pq (while the number of lattice lines 
crossing the side ps is equal to the number crossing the side 
qr) The difference in the number of lattice lines leads to a clo- 
sure failure of the Burgers circuit. This closure failure implies 
that the square contains an excess of dislocations of the same 
sign, the sum of the length of their individual Burgers vector 
is equal to the length of the total Burgers vector. Note that 
the presence of a dislocation is indicated by a T like symbol. 
The vertical bar of the T symbol indicates the side on which 
the extra half plane of atoms is inserted into the lattice, while 
the horizontal top bar indicates the point at which the extra 
half plane of atoms terminates. Thus the Burgers vector is 
parallel to the top bar of the T. 



In constructing a Burgers circuit we take the same num- 
ber of steps in the horizontal and vertical directions, how- 
ever due to the change in the lattice spacing there will 
be a closure failure on the side sr. Assuming the closure 
failure is due to an excess of dislocations of the same sign 
(which are the origin of the extra half planes), then the 
total Burgers vector, due to all the dislocations enclosed, 
is given by 



B^{r) = Ar - L{r)a{r + Ar), 



(16) 



where L(r) is the number of steps taken in the Burgers 
circuit on the side pq, (note that the direction of the 
Burgers vector is perpendicular to the direction in which 
the density is increasing). From Eq. p6)) we have 

B^{r) = Ar - L{r)a{r + dr) 

Ar 
a(r -|- Ar) 

= Ara{r + Ar)[a^'^(r + Ar) - a~^{r)] 
Taylor expansion to first order yields 

da-i(r) 



a{r + Ar) 



L[r) 



B^{r) = Ara{r)[ 



Ar 



dr 
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to get the Burgers vector density we divide through by 
the area of the square and take the hmit Ar — > 0, giving 



hm , T — a{r)- 



Ar^ (Ar)2 



1 d 
2d^ 



dr 



In p(r) 



(17) 



(Generahzing Eq. (I17p . the Burgers vector density in a 
2D plane is given by the vector field Q 



b(r) = a(r)z x Va ^(r), 



(18) 



where z is the unit vector perpendicular to the surface.) 
The density of the Burgers vector is then 



2i?2 



(19) 



where b^{r) is positive upon substituting Eq. ( llOp into 
Eq. (I17p and negative upon substituting Eq. (I12p into 
Eq. 1T7)1 . This then defines a Burgers vector density field 
present throughout the lattice. The total Burgers vector 
within a radius r can be found by integrating Eq. (1191) 
over the area of the disk. 

We have assumed that the changing lattice density is 
due to the insertion of extra half planes into the lat- 
tice. Thus to get the density of dislocations, we divide 
the Burgers vector density by the distance between crys- 
talline rows, h = a\/3/2, which gives for the hard wall 



Mr) 



N 



(20) 



and there is a similar result for the parabolic case [1]. 
Hence, the number of dislocations within a radius r can 
be found by integrating Eq. ([20|) over the area r' < r. 

In addition to dislocations we expect disclinations in 
the lattice. These are present if the Burgers vector den- 
sity is rotational. The density of disclination charge, s(r), 
is given by the z component of the curl of b(r) which re- 
duces in our situation to 



s{r) 



ld_ 

r dr 



irb4,) = -VHnpir)^R{r), 



(21) 



where R{r) is the scalar curvature and is equal to the den- 
sity of disclination charge. This quantity is also known in 
the theory of plasticity as the incompatibility _J2]. This 
relationship can be generalized to include free disclina- 
tions: 



s(r) = s(r) - e,fcVfc&,(r) = R{r) 



(22) 



In analogy to a dielectric the s(r) term is the charge den- 
sity of the free disclinations, induced say as a consequence 
of the topology of the space in which the lattice is em- 
bedded, and — eifcVfc&i(r) is the polarization contribution 
from dislocations 31. It turns out that for the hard wall 



case the six disclinations induced by the disk topology 
are always at the edge of the system. For the harmonic 
potential they are close to the edge in small systems, but 
as N increases, the six topologically induced disclinations 
migrate towards the interior. Thus for values of r away 
from the hard wall, the density of disclination charge, for 
the hard-wall system, can be found by substituting Eq. 
(Uni) into Eq. (EI]), 



1 



-V^ Inp(r) 



1 



1 



i?2 



1 



(23) 



Integrating Eq. ( [23]) over the area r' < r gives the 
total disclination charge within a radius r 



E(r) = 2tt r s{r )dr 
Jq 

/ 7' \ 2 TT 



(24) 
(25) 



If the lattice only contains disclinations with charge ^, 
then the number of disclinations within a given radius is 



Ndisc{r) 



S(r) 
tt/3 



(- 



2 ' 



(26) 



where for the hard wall confined system we expect the 
lattice interior to contain an excess of 7 coordinated 
disclinations. For the harmonically confined system the 
number of 5 coordinated disclinations induced by the 
changing density can be similarly calculated and equals 
Ndisc{r)- The fact that both the density of the Burgers 
vector Eq. (|19p and the disclination charge density due 
to the changing density are equal but opposite for the 
two systems suggests that they are "mirror images" of 
each other. It makes sense then to compare and contrast 
the properties of these two systems. 



C. Smooth part of the energy 

In the following we examine the various terms which 
contribute to the overall energy of the system of charges 
in a hard wall confining potential. A similar study has 
already been made for the system with parabolic confine- 
ment [2]. We believe that the energy will have the form, 
as N becomes large, of a series in decreasing powers of 

iV2 nI N 7V5 

Esmooth ==1^1-5- + ^2—^ + '*3^ + l^i—^ + K5, (27) 
11 ix H ix 

The first coefficient ki — 7r/4; this is calculated next. 
We have obtained from our numerical estimates of the 
ground state energy of systems with varying values of N 
the following estimates of the other coefficients: K2 = 
-1.562033, K3 = 0.975852, K4 = -0.008196 and K5 = 
-0.307608. 
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The first term is the 'electrostatic energy'. It can be 
calculated by approximating the density by its continuum 
limit form, to give 

upon recognizing the second integral as the electrostatic 
potential /i, see Eq. this can be written as 

The next largest term is the 'correlation energy' which 
is the first correction to the continuum limit approxi- 
mation because the charges are discrete. Koulakov and 
Shklovskii ^ suggested that it could be estimated by us- 
ing the local density approximation (LDA) , which states 
that for a large enough cluster, locally the density can 
be assumed to be constant, so the correlation energy of 
a lattice with non-uniform density should be the same 
as the first correction to the electrostatic energy of an 
infinite system of uniform density: 



Ecorr = -7; 2Trrdrp{r)(3^p{r) 
^ Jo 

where the value of (3 depends on the geometric properties 
of the lattice. For a triangular lattice /3 = /?a = 3.921034 
and for a square lattice (3 = f3a = 3.898598 [13]. Visual 
inspection of the ground states of the system (see section 
IV. A. 5) suggest that locally the lattice is triangular, thus 
in calculating the correlation energy we expect that P — 
/3a • If this expectation were correct it would predict that 
K2 = —1.5642939, which is quite close to our numerical 
estimate —1.562033. We discuss in the following sections 
some possible explanations for this small discrepancy and 
also how the remaining k coefficients were obtained. 



IV. NUMERICAL RESULTS 

In this Section we present the results of our numerical 
simulations and compare them with the continuum model 
developed in Section III. 

Simulated annealing experiments were carried out for 
the system with a hard wall boundary, as described in 
Section II, for systems ranging in size from = 2 to 
A^ — 100. Our results either agreed with or gave a 
slightly lower energy than those published by Kong et al 
0. Simulated annealing experiments were also carried 
out for larger systems, A^ = 150, 200, 250, 500, 1000, 
2000 and 5000. To ensure good results a range of an- 
nealing schedules were tried. Upon finding the optimal 
schedule, minimization was repeated as many times as 
possible, thus not only improving the chances of finding 
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FIG. 2: Histograms of the fraction of the total charge of the 
system within a given radius r. Results for clusters containing 
1000, 2000, and 5000 charges are colored red, green and blue 
respectively. The continuum limit result is given by the black 
dashed line. Note that the step- like behavior towards the edge 
indicates that the charge in this region is concentrated into 
a series of concentric shells; the charge arrangements in this 
area are very different from those in the cluster interior. 



a good result but also generating a collection of states 
which could be used for further analysis. 

For the system with parabolic confinement we restrict 
our efforts to large system i.e. N = 1000, 2000 and 5000. 
Numerical experiments have already been carried out for 
small systems by Koulakov and Shklovskii, see and 0]. 
Starting with a random initial configuration of charges, 
a conjugate gradient algorithm was used to minimize the 
energy of the system. In each case this process was re- 
peated 1000 times and the cluster with the lowest energy 
was identified. 

For each system the best result was triangulated; this 
was done by projecting the charges onto a paraboloid 
and using the Delaunay triangulation package Qhull [l^ . 
The results were then displayed using the graphics pack- 
age Geomview [Tsj . An additional routine was used to 
highlight defects in the clusters, points with five near- 
est neighbors were colored red, while those with seven 
and eight nearest neighbors were colored green and blue 
respectively. 



A. Hard Wall System 

1. Distribution of Charges 

In the continuum limit the charge is distributed ac- 
cording to Eq. (|lip ; this quantity can be compared with 
its actual value in a finite sized cluster which we call 
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0.2 0.4 0.6 0.8 1 

r/R 

FIG. 3: Plot of AA'^ for the three largest systems, scaled by 
%/iV versus distance r from the center. AA^ is the difference 
between the number of charges within a given radius for fi- 
nite systems Nfinir) and the continuum result N{r) (which 
is given by Eq. The black dotted line is an attempt to 

fit the data by replacing by — CN^ in Eq. ( fTT|) with 
1.6. 
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FIG. 4: A plot of Nfin{R) against A''^''^ for systems in the 
range N=150 to 5000. To find Nfin{R) we simply counted 
the number of charges in the outermost shell in each system. 



Nfin(r), where for a given radius r, Nfin{r) is the num- 
ber of charges enclosed within that radius. We choose to 
compare the integrated quantity as opposed to the charge 
density itself as this yields a less noisy result. Fig. [5] gives 
the fraction of the total charge enclosed as a function of 
radius, i.e. Nfin{r)/N, for the three largest systems sim- 
ulated (in each case the result with the lowest energy 
is used). Also shown for comparison is the fraction of 
charge enclosed in the continuum limit, i.e. N{r)/N. By 
scaling the charge enclosed in this manner, different sized 
systems can be easily compared. The curves in Fig. [5] 
suggest that with increasing system size the charge dis- 
tribution approaches the continuum result. 

However, there is a systematic difference between the 
charge distribution in the continuum limit and that for 
finite sized systems, which is not just a local effect but 
varies on the scale of the radius R of the system. Its pres- 
ence is revealed on plotting AN{r) — Nfin{r) — N{r), 
for the three largest systems. There is a correction to 
the continuum expression for N(r), which is of order 
\/N. (The continuum expression for N{r) is of order 
N). AN{r)/VN is plotted in Fig. [3]for the three largest 
system sizes which we studied. It shows that there is a 
deficiency of charge in the system interior. Furthermore 
AN falls to zero only at the edge of the system, which 
means that the missing charge is to be found here in the 
form of extra charges, Ns in number, at the surface. The 
leading term for the density is the result in the contin- 
uum limit and is of order N/ i?^ as in Eq. (fTU]) ; Fig. ^ 
shows that there are corrections to it of order /R^. 
We have been unable to obtain any analytic understand- 
ing of this correction term, but we have noticed that it 
can be quite well approximated by replacing N in Eq. 
([ini) by TV - iV^ where TV, = CN^ and C « 1.6. We 
shall estimate the number of charges in the outer shell of 



the system i.e. on the hard wall, and show that in fact 
the "excess" surface charge N^, (which is of order A^a) is 
only a small contribution to it at large N as the number 
of charges in the outer shell increases as A^3. 

The number of charges we would expect to find at the 
edge of the system can be estimated by the following ar- 
gument: let the lattice spacing at the edge of the crystal 
be given by Oe, then the number of charges contained 
within a disk centered on the origin of radius R — is 

N{R - ae) - TV (^1 - ^1-{1-Xf^ , 

where X — a^/R. Expanding the (1 — XY term to first 
order gives 

N{R-ae) mn(i- \/2x) . 

Therefore the number of charges within a distance tte 
from the edge is 



N -NiR-a^) 



Imposing the condition that A'e must be equal to the 
number of charges on the perimeter Np = 27r i?/ae, yields 
the ratio 



R 



A^ 



(28) 



hence the number of charges on the edge, in the con- 
tinuum limit, scales with Na. To compare with the re- 
sults of our numerical simulations we therefore plotted 
NfiniR), the total number of charges in the outermost 
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FIG. 5: Test of the LDA approximation. If the approximation 

would be independent of Ti/R 
-3.921034. 



were perfect /\E(ri) / ^/pin 
and equal to the numerical value 



shell of each cluster, against N'^^^. As shown in Fig. 2] 
there is good agreement between these estimates and the 
results of our numerical simulations. 



2. The Local Density Approximation 

In Section III.C we calculated the correlation energy of 
the lattice using the local density approximation (LDA) 
which is based on the idea that locally the lattice appears 
to be triangular. Physically this seems to be a very rea- 
sonable approximation, but we observed that the numer- 
ical data was not in perfect agreement with this approx- 
imation. In this Section we shall investigate the matter 
further. 

The correlation energy of a given charge will be de- 
fined as the difference in energy between the interaction 
of a charge with all the other charges in the system after 
the subtraction of the interaction of the charge with the 
continuum approximation for the other charges (Tsj : 



AEin) 



^ 1 

X! Iv... _ 



N 

E 



1 t:N 



v,-vA 2R 



(29) 



In the LDA, this quantity would be expected to equal 
—Pa \fpiTi)- In Fig. ([5]) we have plotted AE{ri)/ \/p{n) 
against ri/R to test this expectation. The agreement is 
non-existent. 

The origin of this discrepancy is not hard to find. It 
arises because AE{ri) is a quantity of size N"^ /R. \fp{j\) 
is of this magnitude, but so are also the contribution of 
the deviations of the density from p^{r) which are of 
order /R^ . Alas, we do not have a calculation of these 
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FIG. 6: Histogram giving the total number of disclinations 
within a given radius. Clusters containing 1000, 2000, and 
5000 charges are colored red, green and blue respectively. The 
analytical curve is given by the black dashed line. The inset 
is a magnification giving the last 20% of the graph showing 
there is good agreement right up to the edge. 



deviations. However, we have found that our attempts to 
model them by changing N in the expression for (r) to 
N — CyJ\N) and allowing for the compensating surface 
charge needed for charge neutrality does at least reduce 
the discrepancy in Fig. So it could be that with 

proper allowance for these corrections to the density, the 
LDA might still be valid. 



3. Distribution of Disclinations 

The number of excess disclinations located within a 
given radius is predicted by Eq. (|26| . To compare this 
with our simulations we use the following method: for a 
given radius, we count the number of positive and neg- 
ative disclinations enclosed, where we expect from Eq. 
((26)) that in the interior of the lattice the number of neg- 
ative disclinations will always be greater than the number 
of positive disclinations. We adopt the convention that a 
7 coordinated point counts as one negative disclination, 
an 8 coordinated point as two negative disclinations and 
a 5 coordinated point counts as one positive disclination. 
In Fig. [5] we plot the number of excess negative disclina- 
tions within a given radius and compare with Eq. (I26p . 
It is evident that there is convergence with increasing 
system size. For every negative disclination in the lat- 
tice interior there is a compensating positive disclination 
on the lattice edge. In addition, the lattice edge also 
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FIG. 7: ^ against iV'i plotted for N=150 to N=5000. The 
red line is a fit of the form y=0.975852x + 0.00223295 




contains 6 extra positive disclinations due to Euler's the- 
orem. 



FIG. 8: A plot showing the difference between the energy for 
small clusters found by simulations and the smooth part of 
the energy. 



4- Cluster Energies 



In this section we compare the scaling of the energy as 
given in section III.C with the energies found by simu- 
lated annealing. 

Having identified the electrostatic and correlation en- 
ergies as the first two leading order contributions to the 
energy of the cluster, which scale as N'^ and N^-^ respec- 
tively, we make the ansatz that the remaining terms in 
the series descend in powers of ^/N, thus 



Es 



mooth 



N 



Ki—fT- + K5- 
it 



(30) 



The second term is divided into two parts: /3a is the 
correlation energy if the lattice is strictly triangular ev- 
erywhere and the LDA applies; the term containing the 
parameter a is a correction to this assumption. 

We now compare the scaling of the energy given by 
Eq. ( 130]) with the energies found by simulated anneal- 
ing. Since the electrostatic and correlation/surface terms 
dominate for large systems, it is natural to ask what in- 
fluence the remaining terms may have, so we plot 

AE if TTiV^ n~iVi\ , , 

against N~i, see Fig. [71 From the intercept we find 
a = 0.00223295 and from the gradient K3 = 0.975852. 
Let us replace the coefficients multiplying the N^-^/R 
term in Eq. pOp with 



K2 



(3o 



Pa, 



(32) 



where /3o — 3.915436 is the 'observed' value of /5. It dif- 
fers from /3a by approximately 0.15%. This might be due 
to a failure of the LDA approximation, but it might also 
be due to the fact that we are almost certainly not ob- 
taining the ground states for each value of N, which has 
uncertain consequences for the accurate determination of 
the coefficient K2. 

There is no appreciable bending of the curve shown 
in Fig. ([71), which suggests that the remaining terms in 
Eq. ( I30p are insignificant in the large N limit. (The 
coefficients K3, K4 and K5 were actually determined by 
using a curve fitting package in the range iV = 2 to 100). 
Furthermore as all the data points lie almost perfectly 
on a straight line, we conclude that the numerical results 
are fairly reliable. In fact for large systems, this plot 
often enabled us to determine which clusters generated 
by the simulated annealing algorithm required further 
annealing. These were the clusters for which the data 
points were slightly above the fit plotted in Fig. [71 

Having determined the form of the "smooth" part of 
the energy we can determine the fluctuating part of the 
energy, which is defined as 5E = E — E smooth- The 
results are shown in Fig. [51 Some clusters have a par- 
ticularly low value of 5E such as = 34, 49 and 62. 
These clusters are known as "magic number states" and 
posses a high degree of geometrical symmetry. To illus- 
trate this point Fig. ^ shows a Delaunay triangulation 
of the system with 62 and 69 charges, which correspond 
to the lowest and highest points, respectively, in Fig. ([5]). 



10 




500 Charges 1000 Charges 



FIG. 9: Clusters with a hard wall confining potential. The cluster with A'^ = 62 has the lowest measured value of SE, 
while the cluster with A'^ = 69 has a particularly high value. This is because the N — 62 cluster possesses a high degree of 
circular symmetry, which is commensurate with the hard wall boundary. The lack of circular symmetry in the A'^ = 69 cluster 
is evident from the off-center negative disclination. In addition the boundary imposes a shell-like order towards the edge of the 
cluster, which is indicated by the presence of a large number of seven and eight coordinated disclinations. For larger systems 
the shell-like order is confined to within a few lattice layers of the edge, the lattice interior is composed of a triangular lattice 
with some charged grain boundaries. 
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5000 Charges 

FIG. 10: Clusters with a hard wall confining potential. The remarkable feature of this system is the arched like structure 
towards the edge. This curvature is due to an excess of negative disclinations in the lattice interior. Most of the charged grain 
boundaries appear to be located in the region between the relatively undistorted triangular lattice and the arches. We believe 
that the grain boundaries act to screen off the two incommensurate regions. 



5. Small, Medium and Large Clusters 

Small clusters, those which contain less than 100 par- 
ticles, are dominated by the circular hard wall boundary. 
In particular, clusters containing less than 50 charges 
show no hint of hexagonal order. Instead the charges are 
arranged in concentric shells, where many of the charges 
have 7 or 8 nearest neighbors. We identify this type of 
ordering as "shell- like" ; it is present at the edge of all the 
systems. For systems containing more than 50 particles, 
a region of the cluster free from the influence of the hard 
wall boundary begins to emerge; this region looks like a 
triangular lattice. At this point there is a competition 
between shell-like order and the emerging hexagonal or- 
der. The clusters with the lowest energy are the ones 
which can satisfy both simultaneously. Both A'' = 62 and 

= 69 are clusters which contain a region where most 
of the charges have six nearest neighbors. The charges 
near the center in both clusters can be considered to be 



part of a disclination as the central charge has a different 
coordination number compared to its neighbors. Of all 
the clusters, A^ = 62 has the lowest value of 6E while 
A^ = 69 has the highest. The difference is that for N=69 
the disclination is off-center which breaks circular sym- 
metry while in the case of A" = 62, the disclination is 
perfectly in the center of the cluster and is surrounded 
by a ring of charges which have six nearest neighbors. 

For the medium sized clusters, i.e. 100 < A^ < 1000, 
the influence of the hard wall boundary on the clus- 
ter is contained within a narrow annular region at the 
edge. This again is the shell-like region and contains 
a large number of seven and eight coordinated disclina- 
tions, which have the effect of destroying the hexagonal 
order. On the other hand, the internal region of the lat- 
tice is comprised of a relatively undistorted triangular 
lattice lattice. The lattice interior contains an excess of 
negative disclinations. However, for every excess negative 
disclination in the interior there is a compensating posi- 
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r/R 

FIG. 11: Plot of AN for the three largest systems with a 
harmonic confining potential, scaled by JV°'^ versus distance 
r from the center. 



live one on the lattice edge. Using the rule that a charge 
with 7 or 8 nearest neighbors is a disclination with topo- 
logical charge -1 and -2 respectively, while a charge with 
5 nearest neighbors has charge +1, the total topological 
charge of all the clusters was found to be +6. The 6 posi- 
tive disclinations due to Euler's theorem are also located 
on the lattice edge, hence Euler's theorem is satisfied. 

Also with increasing system size isolated disclinations 
become less common and are replaced by small topologi- 
cally charged grain boundaries; these are chains of alter- 
nating positive and negative disclinations, which contain 
in total one excess negative disclination. As shown in 
Fig. , the alternating Burgers vectors of the charged 
grain boundary cancel out and the overall arrangement 
constitutes a disclination. 

Large clusters are those which contain more than 1000 
charges, see Fig. PH)) . In this regime, the interior con- 
tains large areas of lattice separated by charged grain 
boundaries, which are numerous and long. As will be 
shown in Section V.A, excess disclinations of one sign in 
the lattice interior generate lattice curvature, the effect 
of which is to bend the lattice lines into a series of arches. 
This can be seen with great clarity in the N = 5000 sys- 
tem. For this system, the bending is sufficiently strong 
to make the arched structure incompatible with the or- 
dinary triangular lattice. It would seem that the grain 
boundaries arrange themselves to screen off the two re- 
gions. 



B. Harmonically Confined System 

In this Section we consider a 2D cluster of charges con- 
fined by a harmonic confining potential. This system 
has been studied in considerable detail by Koulakov and 




FIG. 12: A plot of the number of excess positive disclina- 
tions within a distance r/R where we use the rule that a five 
coordinated and a seven coordinated disclination count as a 
single positive or negative disclinations respectively. Results 
for clusters containing 1000, 2000 and 5000 clusters are col- 
ored red, green and blue respectively. The analytical curve 
is given by the black dashed line. Note that at the edge of 
the cluster the number of excess positive disclinations falls to 
+6. Unlike the number of excess disclination in the hard wall 
system, the number in this system lags behind the continuum 
value significantly. 



Shklovskii, see 

m and a, and here we expand on their 

work. 

As for the system with the hard wall confining poten- 
tial, the actual density differs from the continuum limit 
result by subdominant contributions which are difficult 
to quantify. Plotting AN = A^/m(r) - N^ir) for the 
three largest systems shows that there is a deficiency of 
charge in the system interior which is compensated for 
by a shell of excess charge at the edge of the system. A 
good collapse of the data for different values of N was 
found as in when we plotted 

AiV ^ iV/.„(r)-iV^(r) 

NO-e ]\[0.6 ' ^'^'^) 

but the origin of this correction with the power 0.6 is a 
mystery to us. 

In the continuum limit, the density of excess disclina- 
tions and the density of the Burgers vector for this and 
the hard wall systems are the same except for the sign. 
It follows that the number of excess disclinations within 
a given radius is the same for the two systems. A com- 
parison with our numerical simulations is given in Fig. 

Koulakov and Shklovskii demonstrated that provided 
the cluster is small {N < 150) then there exist cer- 
tain magic number states which contain only 6 five- 
coordinated disclinations demanded by Euler's theorem. 
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FIG. 13: Clusters with a harmonic confining potential. 1000 charges: This cluster contains both charged grain boundaries 
and also some isolated dislocations. These dislocations are not involved in screening, but serve to relax the strain energy. 2000 
charges: With increasing system size, the charged grain boundaries continue to grow in length. Also the isolated dislocations 
become more numerous, especially near the edge of the system. Note these dislocations are oriented with the 5 coordinated 
disclination pointing towards the center of the system. There is a rapid rise in the number of positive disclinations just before 
the lattice edge which are canceled by the negative disclinations close to the edge. 5000 charges: For very large systems, the 
additional dislocations required to accommodate the changing density form uncharged grain boundaries (i.e grain boundaries 
with no over-all topological charge). Clear examples of such grain boundaries can be seen towards the edges of the cluster. 



14 



These disclinations are not right on the lattice edge but 
are close to it. Beyond this limit {N > 150) the discli- 
nations are always accompanied by a screening cloud of 
dislocations. It was further shown that provided the clus- 
ter is not too large {N < 700) then the only defects 
in the lattice are the disclinations and their screening 
clouds, which together form 6 small separate topolog- 
ically charged grain boundaries. In fact provide that 
N < 700, the lattice can be divide into an inner region 
and an outer region, where the boundary between the 
two is given by the radius at which the charged grain 
boundaries are located. The inner region is a largely un- 
deformed triangular lattice while in the outer region the 
lattice lines are curved. 

Beyond this limit {N > 700) dislocations not associ- 
ated with screening begin to appear, these dislocations 
are present in order to ensure that Eq. (f23|) is satisfied, 
consequentially the lattice cannot be divided so neatly 
into two separate regions [3]. 

The numerical work of Koulakov and Shklovskii only 
dealt with systems containing less than 700 charges. Our 
work is concerned with the behavior of the system as 
the cluster grows beyond this limit. An examination 
of the cluster containing 1000 charges, see Fig. p^ . 
shows that in addition to the topologically charged grain 
boundaries the system also contains a few isolated dislo- 
cations. These dislocations are orientated so that the 5 
coordinated disclination points towards the center of the 
lattice while the 7 coordinated disclination points radi- 
ally outwards. With increasing system size the general 
trend is that the screening cloud around the disclinations 
continue to grow in length. In addition the isolated dis- 
locations become increasingly numerous, see the cluster 
containing 2000 charges in Fig. p^ . Eventually in ad- 
dition to isolated dislocations we can observe uncharged 
grain boundaries (i.e extended chains consisting of alter- 
nating disclinations but which have no overall topologi- 
cal charge), see the system containing 5000 charges. Fig. 

Even though the continuum disclination charge density 
is the same for the hard wall system and the system with 
harmonic confinements, nevertheless there is a remark- 
able difference between the two in their approach to the 
continuum limit with increasing system size. Compared 
to the hard wall system, the number of disclinations in 
this system is far fewer. From Fig. (fT^ we can see that 
the number of disclination within a given radius lags be- 
hind the continuum value, while the opposite is true in 
the hard wall system. For instance, for 5000 charges in a 
harmonic trap the number of excess disclinations at the 
edge is about 25, while in the hard wall case this is 400. 
Even though the curvature of lattice lines in the con- 
tinuum limit (which is given by the density of Burgers 
vector) is the same for the two systems, the curvature of 
the lattice lines appears to be far less pronounced for the 
system with harmonic confinement. This may be due to 
the lack of disclinations, since the curvature depends on 
the number of disclinations enclosed within a radius, see 




FIG. 14: A close up of a small grain boundary with a to- 
tal topological charge of -1. The arrangement can also be 
thought of as a pair of dislocations with opposite Burgers 
vectors which cancel. Drawing a Burgers circuit around the 
grain boundary leads to a closure failure which, like a discli- 
nation of charge -1, increases with distance from the center. 



Eq. HMD. 

Unlike the hard wall system, which requires an ex- 
cess of 7 coordinated disclinations in the interior, the 
harmonic system already contains six free 5 coordinated 
disclination due to Euler's theorem. In the case of the 
hard wall system these topologically induced disclina- 
tions are pushed to the edge of the system while in the 
harmonic case these disclinations sink into interior and 
help in accommodating the non-uniform density. In con- 
trast to the hard wall system, the density of free discli- 
nations, given by s(r) in Eq. f l22p . cannot be ignored. 
From Fig. it can be seen that up to r = 0.8R the 

total disclination charge does not exceed -1-6. Thus up 
to this point the mechanism by which the lattice adapts 
to the decreasing density towards the edge of the system 
depends on the arrangement of these free disclinations. 
In order for the total disclination charge to match that 
given by Eq. dislocations arise to screen the discli- 

nations. These screening dislocations have the effect of 
smearing out the disclination charge. To explain what 
we mean by this, consider the following qualitative argu- 
ment. To define a disclination we must be able to draw 
a Burgers circuit around it. The smallest Burgers circuit 
around a single disclination has a radius equal to the lo- 
cal lattice spacing. This then defines the minimum size 
of the defect. By screening the disclination to form a 
charged grain boundary, of the type shown in Fig. (I14[) . 
the radius of the Burgers circuit needed to enclose the 
object becomes larger, hence the size of the defect is in- 
creased. However, the total disclination charge enclosed 
is still the same. Thus the total density of disclination 
charge, i.e. s(r), is reduced. Or to put it another way, 
unlike an isolated disclination, which is a single point in 
the lattice, the charged grain boundary is an extended 
object. However, the total charge of the grain bound- 
ary is still ±1. Thus we can consider this charge to be 



15 



smeared out over its length. 

Beyond r = 0.8R there is a sharp increase in the num- 
ber of dislocations, which for the largest systems con- 
dense into uncharged grain boundaries. These disloca- 
tions are oriented such that the 5 coordinated disclina- 
tion points inwards while the 7 coordinated disclination 
is on the lattice edge. Thus towards the edge of the 
lattice there is a sudden jump in the number of excess 
positive disclinations followed by an equally sudden fall. 
Perhaps as the number of charges in the system is in- 
creased further, these dislocations might unbind so that 
the 7-coordinated disclinations are pushed to the lattice 
edge while the 5-coordinated disclinations remain in the 
interior. 



V. DISCUSSION 
A. Lattice Curvature 

The most remarkable feature of the numerical simu- 
lations is the bending of lattice lines towards the edge 
of the cluster; which is particularly strong in the larger 
systems - see for example Fig. (|10|) . This bending is due 
the fact that the lattice contains an excess of disloca- 
tions of one sign. The original argument explaining this 
phenomena was given by Nye p^ . who showed that the 
curvature of the lattice is equal to the Burgers vector den- 
sity. However, Nye's exposition assumes that the lattice 
spacing remains constant everywhere. Blindly applying 
Nye's result to our system would predict that the lattice 
lines bend in the opposite direction. With a slight mod- 
ification we can adapt Nye's argument to explain lattice 
curvature in crystals with a changing density. 

Consider a square section of lattice pqrs of dimension 
Ar as shown in Fig. P^A). the lattice has a constant 
density everywhere and so there are an equal number of 
lattice lines crossing each side of the square. If on the 
other hand the density is increasing in the r direction it 
means there must be more lattice lines crossing the side 
sr than pq. Drawing the Burgers construction around 
the square leads to a closure failure. Consequentially 
the square must contain an excess of edge dislocations 
of one sign, the sum of their individual Burgers vector 
being equal to the total Burgers vector B. This situation 
is shown is Fig. (fTSl B). it can be imagined that the total 
Burgers vector is split into two equal parts and all the 
dislocations are contained within the triangles p't's' and 
q'u'r'- both of which subtend an angle of 6'/2 - it follows 
that 



(34) 



where for small values of 6 we can ignore higher order 
terms. It is useful to imagine that the original lattice has 
been plastically deformed into the quadrilateral p'q's'r'. 
However, this picture is misleading as none of the lattice 
lines which were originally parallel to ps and qr suffer a 



change in length; thus the state shown in Fig. (fTSlB) 
is for illustrative purposes only. The true state of the 
deformed lattice is shown in Fig. (fTSlC). By mapping 
the sides pq and sr of the original square onto the circular 
arcs p 'q ' and s 'r ', all of the lattice lines originally parallel 
to ps and qr remain of length Ar - at the expense of being 
no longer parallel to each other. To find the curvature k 
of the bending of the arcs, let the length of the lines o'p' 
and o'q' be equal to L. We use the relationship 



Ar 



where k — 1/L. Upon substituting for 6 from Eq. 
and taking the continuum limit we have 



k(r) 



lim 

Ar^ 



nB 



b(r). 



(35) 



Hence the curvature is equal to the Burgers vector den- 
sity. If on the other hand the lattice density is decreasing 
in the r direction then we expect the sense of curvature 
to be reversed. Furthermore, for lattice lines not paral- 
lel to the local Burgers vector, then the curvature of the 
lattice line ko depends on the angle 7 it makes with the 
local Burgers vector [iBl 



ko(r) — k(r) cos 7. 



(36) 



The next logical step is to show that the lattice cur- 
vature in the systems which we have simulated is given 
by the Burgers vector density. Consider the deformed 
hexagon shown in Fig. (|16p where the direction of b is 
marked by an arrow. We expect the curvature of the 
lattice lines to be described by Eq. f [36|) . To make a 
connection with our simulations, for the N = 5000 sys- 
tem with a hard wall boundary shown in Fig. (jlOp. each 
hexagon is decomposed into three arcs, each of which can 
be further decomposed into three points as in Fig. (fTBjl . 
By fitting the points to the equation of a circle we deter- 
mined the curvature of each arc. We assume that each 
arc yields the curvature of the lattice at the center of the 
hexagonal cell. Depending on the radial distance of the 
cell from the center of the lattice we expect this curva- 
ture to have any value between and |b(r)| (depending 
on the orientation of the lattice line with respect to the 
local Burgers vector density field) . For the hard wall sys- 
tems with 5000 charges. Fig. (IT7)) shows the curvature 
of each such arc against radial distance; also plotted for 
comparison is the Burgers vector density given by Eq. 
([TO]) . Ideally the curvature ought not to exceed the limit 
set by Eq. (I19p. however, this is not possible in reality 
as lattice lines close to disclinations suffer much greater 
curvatures than Eq. (|19p would allow. From the general 
similarity of the curvature data and Eq. I19p we conclude 
that the cause of the bending of lattice lines is indeed due 
the plastic deformation of the lattice. 

Next we wish to make the connection between curva- 
ture as defined by the work of Nye to some well-known 
results about parallel transport of a unit vector about a 
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FIG. 15: [A]: An undeformed lattice of dimensions Ar. [B]: A lattice of increasing density, where the change in lattice density 
depends only on the r direction. In effect the original square lattice has been plastically deformed into the quadrilateral p'q'r't'. 
This diagram is geometrically incorrect as only the lines originally parallel to pq and sr should suffer a change of length, while 
those originally parallel to ps and qr should remain of the same length. [C]: The true deformed state of a lattice with increasing 
density, the lattice lines pq and sr are deformed into the circular arcs p'q' and s'r'. This state is geometrically correct as all of 
the lattice lines which were originally parallel to ps and qr are still of the same length. Note this construction only tells us how 
to calculate the curvature of lattice lines which were originally parallel to B 




FIG. 16: This diagram shows a hexagonal cell in a lattice with 
a non-zero Burgers vector density field b. The hexagonal cell 
can be decomposed into three arcs which cross the center of 
the cell, if the arc makes an angle 7 with the vector b then we 
expect its curvature to be given by |b|cos7. Thus as shown 
in this diagram the horizontal arc which is perpendicular to 
b has no curvature. 



vector in terms of the disclination charge enclosed 
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(37) 

If as before we let E(r) be the total disclination charge 
within a disk of radius r, then using Eq. ( I35p we can 
write Eq. ((37|) as 



27rr 



(38) 



This equation expresses the fact that the curvature at a 
distance r from the center of the system depends on the 
total disclination charge enclosed within a disk of radius 
r. For a given curve the curvature is defined as the rate 
of change of the angle of its tangent vector ui, thus for 
the circular path enclosing a disclination charge S(r) we 
have 



1 duj{r) 



disclination \1T\ . If a unit vector is transported on a lat- 
tice along some path enclosing disclinations, then when 
the vector returns to its original position its orientation 
will have changed. The total change depends on the num- 
ber of disclinations enclosed and in a triangular lattice 
must be a multiple of 7r/3. To show this connection we 
can invert Eq. ipT]) to express the density of the Burgers 



To find the total change in the angle of the tangent vector 
we integrate over the length of the circular path giving. 
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27r 



rk^{r)d4> = i;(r). 



where in a real lattice the disclination charge is quan- 
tized. This result explains why the orientation of the 
lattice cells in images such as Fig. ([TU|) is observed to 
rotate upon traversing a circular path centered on the 
origin. 

It should be noted that as long ago as 1955 the geom- 
etry of imperfect lattices had been developed by Kondo 
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FIG. 17: Curvature of lattice lines in the system with 5000 
charges in a hard wall confining potential. As discussed in 
the text, each hexagonal cell can be decomposed into a set 
of 3 arcs; for each such arc the curvature is computed and 
plotted against the radial distance of the parent hexagon - 
this is shown by the black dots. Theoretically the maximum 
curvature that an arc can have is given by |b(r)| which is 
shown by the red dotted line. Thus depending on the orien- 
tation of the lattice line, the curvature ought to range from 
to |b(r)|. The fact that the lattice curvature exceeds |b(r)| 
for some arcs is due to the fact that close to the core of a 
defect the curvature becomes much larger than that set by 
the continuum limit calculation. 



and co-workers into non-Riemann differential geometry 
[Tst - By analogy to general relativity, one can think of 
an undeformed lattice as a region of fiat space, which 
becomes warped in the presence of defects. The authors 
demonstrated that dislocations and disclinations gener- 
ate torsion and curvature respectively. 



B. Conformal Crystals 

The systems we have discussed thus far have been 
constructed using optimization algorithms to generate 
ground state configurations. Surprisingly there exists an 
unusual class of 2D lattices which have a non-uniform 
density but which can be constructed by a purely ana- 
lytical method. These structures are known as conformal 
crystals [l^ . As the name suggests the positions of the 
lattice sites can be obtained in the image plane by ap- 
plying a conformal transformation to a regular lattice in 
the z-plane. 

An example of a conformal lattice with circular sym- 
metry is shown in Fig. (|18p. As evident from the way 
in which the lattice lines curve towards the edge and the 
increasing density, there is a resemblance between the 
conformal crystal and the clusters studied in Section IV. 
Initially it was hoped that these conformal crystals might 
help explain the origin of the lattice curvature observed 
in the other 2D systems, but this did not turn out to be 
the case. In this section we show that conformal crystal 
can be regarded as a giant disclination. 



Originally the idea of conformal crystals was used to 
describe a structure formed by a cluster of mutually re- 
pelling magnetized spheres dubbed "gravity's rainbow" . 
In an experiment, metallic spheres were confined to a thin 
rectangular box which is placed in a magnetic field, the 
field induces a magnetic moment in each of the spheres 
which causes them to repel. Under the action of grav- 
ity the spheres crystallize into a lattice with non-uniform 
density, which consisted of a series of arch-like structures. 
The authors suggested that the unusual lattice could be 
obtained by a conformal transformation of a regular tri- 
angular lattice (20| . 

Consider for example a regular triangular lattice in 
the z-plane such as that shown in Fig. ([18]- left). By 
applying an analytical transformation w = f(z) the cor- 
responding coordinates in the w-plane are 

u + iv — f{x + iy), 

where w = u + iv = re'^ and z = x + iy = r'e**^. In 
the case of f{z) ~ z^ the result of the transformation is 
shown in Fig. p^ . This transformation belongs to a set 
of transformations 



= C 



or w 



(39) 



which yield conformal lattices of circular symmetry. Here 
we shall only concern ourselves with the first of these 
transformations. 

Conformal transformations have three important fea- 
tures. Firstly they are locally angle preserving (isogonal) 
transformations. This means that upon mapping an in- 
finitesimal hexagonal lattice cell in the z plane to the 
image plane, the cell will still have all its internal angles 
equal to 7r/3. It is important to note that because of the 
local nature of the angle preserving property this is only 
strictly true for an infinitesimal lattice cell. For a real 
lattice with a well-defined lattice spacing, the mapping 
distorts the shape of the hexagonal cells. This distor- 
tion is stronger towards the center of the lattice than the 
edge, see Fig. (IT51) for an illustration of this. At the origin 
the conformality of the transformation breaks down com- 
pletely. In the case of the transformation w{z) = z^^^, 
instead of preserving angles, the angle of tt/S at the ori- 
gin is halved in the image plane. No matter how small 
the lattice cell enclosing the center in the z-plane is made, 
this breakdown of conformality will remain. If the confor- 
mality of an otherwise conformal mapping breaks down 
at a particular point, then that point is called a critical 
point of the mapping. The critical points of a conformal 
transformation exist at any point where either \dw/dz \ or 
its inverse is equal to zero [2l|. Secondly, a lattice with 
constant density pz in the z-plane will have a density [Toj 



dw 



dz 



(40) 



in the w-plane. Thirdly, it was shown that if the transfor- 
mation is conformal then the density of lattice points in 
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FIG. 18: Left: a circular cut out of a triangular lattice in the z plane, where the angle between each successive red line is 
7r/3. Right: the effect of the transformation w = z2 in the image plane, where the lattice has been produced by allowing 
the range of the z plane to extend to < 4tv. Conformal transformations are locally angle preserving transformations, 
therefore an infinitesimal hexagonal lattice cell in the z plane is not deformed upon mapping to the image plane. However, 
this angle preserving property of the transformation breaks down at certain points which are known as critical points. For the 
transformation w = Z2 this occurs at the origin, where instead of preserving the angle n/3, it is halved. Note that the lattice 
cells are more distorted closer towards the critical point. 




FIG. 19: A coifformal lattice generated by the transformation 
w = 2®''^ which is similar to a disclination with charge -1. Also 
shown are a series of Burgers circuits enclosing the central 
point. 

the image plane is constrained by the following condition 

M 

V2lnp(r) = 0, (41) 

which interestingly is also the condition that the lattice 
has no disclination charge induced by a varying density, 
see Eq. (1^ . 

Thus a circularly symmetric conformal crystal such as 
the one shown in Fig. (|18p is locally hexagonal but has 



curvature. It has an increasing density and except for the 
point at the center has no apparent internal defects (by 
which we mean that a triangulation only shows internal 
points which have six nearest neighbors) . Yet the frame- 
work developed in Sections 2.3 and 2.4 demonstrates that 
a change in the lattice density must be accompanied by 
lattice defects, which in turn are responsible for curva- 
ture. How are these two seemingly conflicting statements 
to be resolved? A clue is provide by the point at the cen- 
ter of the lattice with the anomalous coordination num- 
ber. 

Consider the conformal lattice generated by the trans- 
formation 

w{z) — Z7 , 

which is shown in Fig. (|19p . Superimposed on top of the 
conformal lattice are a number of Burgers circuits which 
enclose the central point. Whereas previously we were 
dealing with a continuum, in which the defects were as- 
sumed to form a gas throughout the system, here there 
is only a single defect at the center of the system. Since 
there are no other defects, any Burgers circuit which does 
not enclose this central point will close. As shown the cir- 
cuits start at a,b... and end at a',b'..., in each case there 
is a closure failure, which increases with distance from 
the central point; indeed the successive Burgers circuits 
trace out a wedge. Thus we propose the central point in 
a conformal lattice can be thought of as a "disclination" . 
The same rules which apply to a disclination also apply 
here, i.e. the wedge angle is quantized by the symmetry 
of the lattice and the wedge in turn can be decomposed 
into a series of half planes [llj . To see how a wedge can 
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be decomposed into a series of half planes, note that the 
difference in the closure failure between any two succes- 
sive Burgers circuits, such as aa' and bb', is always one 
lattice spacing, which implies that between a pair of cir- 
cuits an extra half plane has been inserted at the points 
labeled a,b... The question of where the half planes have 
been inserted is arbitrary as it depends on where the 
Burgers circuits start and end. Thus as for a disclina- 
tion the relationship given by Eq. (|14p also holds for a 
conformal lattice. In the case of = (also for a -1 
disclination) the wedge angle is given by f2 = 7r/3. By 
considering a series of circular concentric circuits whose 
origin coincides with the critical point of the transforma- 
tion, it is possible to draw a series of Burgers circuits and 
calculate the Burgers vector density using the approach 
outlined in section III.B, in either case the result is still 
given by Eq. P7|) . 

Since at the critical point the lattice ceases to be con- 
formal, Eq. I|4ip is true everywhere except for the origin. 
This suggests that it can be written as 

V^lnpM = (42) 

where v is an undetermined constant and is the two 
dimensional delta function. To find we are going to 
assume that the point at the center of the lattice is a 
disclination with charge s(r). We have the following re- 
lationship between Eq. (jH]) and Eq. (dH) 

2S(r) = Inp(r) = i^(5(2)(r). (43) 

Substituting w{z) — Czx into Eq. (HD|) gives 

P{r) = PzX^r^^ = , (44) 

where we have used the relationship r — r x . Thus we 
can write Eq. as 

2S{r) =2{x~l)V^\nr ^ i^S'-^\r), (45) 

Recognizing V^lnr = 27r(5^^^(r) as the two dimensional 
Green's function and canceling out the factor of 2, Eq. 
becomes 

thus the central point of the conformal transformation 
has a disclination charge density of i^' = i^/2 = 27r(x — 1). 
The total disclination charge contained within the disk 
can be found by simply integrating over its area, thus 

E(r) = / / s{r)rdrd9 
Jo Jo 

= 27r(x- 1) / / d^^Hr)rdrd0 
Jo Jo 

= Mx-l), (46) 

In the case of x = 2 which corresponds to the transforma- 
tion w{z) = z^/^, the central point of the transformation 
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FIG. 20: Curvature of the lattice lines for a conformal crystal 
generated by the transformation ■w{z) = 2:2. The red line 
gives the maximum curvature using Eq. I |47|l . The black dots 
give the actual curvature of the lattice line, as calculated by 
the method outlined in Section 4.5. 



has a disclination charge of 6{tt/3) — 2tt, i.e. a discli- 
nation consisting of 6 wedges, each of which subtend an 
angle of tt/S in the z plane. 

This realization that the central point in a conformal 
lattice is actually a disclination can be used to calculate 
the lattice curvature. The conformal lattice shown in Fig. 
(|18p was generated by applying the transformation w — 
02 to a hexagonal lattice, the original lattice contained 
Nz points within a disk of radius R, thus pz = Nz/ttR'^, 
setting X = 2 in Eq. (|44|) yields 

Pw{r) = -ir^Pz, 

Using the framework developed in Section 2.4 we know 
that the maximum curvature k is related to the density 
of the Burgers vector by Eq. ((35|) . thus 

Hr) = |b(r)| = ^|:lnp.(r) = ~^^^r'pz = \, (47) 

where we have used Eq. ()17|) . On the other hand, using 
Eq. ([13]) with = 27r gives |B(r)| (27r)r. Inverting 
the relationship given by Eq. (PT|) yields 

Hr) = |b(r)| = ^^\^{r)\ = i^^'^^r = - (48) 
ZTrr dr Zirr dr r 

Alternatively, using Eq. (I38p . we can assume that the 
disclination charge enclosed is equal to 2t: and this gives 
the same result. Thus a comparison can be made between 
this analytical result and the actual measured lattice cur- 
vature (just like we did for N = 5000 in the hard wall 
case). Fig. pO|. There is perfect agreement between the 
two curves which leads to the conclusion that a conformal 
lattice can be thought of as a type of disclination. 

C. Experimental Realizations 

The system with the hard wall confinement could po- 
tentially be realized in a number of different physical 
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contexts. One possible experimental situation involves 
a collection of polystyrene beads on a disk, located on 
the edge of which is an insulated boundary. The beads 
are then exposed to ionizing radiation and the liberated 
electrons are sucked out of the system by an electrode, 
leaving behind a collection of positively charged beads. 
The cluster can then be annealed to the ground state by 
shaking the disk. 

Alternatively the system could be realized using col- 
loidal particles. For these systems one would have to en- 
sure that the Yukawa screening length was made larger 
than the system size. By using less polar organic solvents, 
screening lengths as long as 12/im can be achieved [2^ . 
Thus to realize large systems one would need colloidal 
particles with diameters of, say, 0.1/im. 

The system with parabolic confinement is of consid- 
erable interest to the field of quantum dots. For review 
articles discussing the fabrication of quantum dots and 
the harmonic confining approximation, see p3| and (23 | 
respectively and the references contained. 

D. Other Systems 

There are a number of systems in which lattice curva- 
ture is quite prominent; one such system is the growth 
of crystals in amorphous films [25l |. Up until now the 
suggestion has been that these systems possess a confor- 



mal geometry because the lattice lines are curved. This 
work demonstrates that this is not necessarily so. The 
condition that a system has a conformal geometry is very 
strict, i.e. the lattice density has to obey Eq. ( I44p . It 
is possible that these systems contain an excess of discli- 
nation charge in the interior, which in turn results in 
the bending of lattice lines. It would be interesting to re- 
examine such systems in light of the results of this paper. 
Even the original experiment which sparked the inter- 
est in conformal crystals, the so-called gravity's rainbow 
structure, has been shown not to be stable [26|, meaning 
that it is likely that the system does not form a perfect 
conformal crystal. Numerical simulations suggest that 
the system is composed of domains which are separated 
by defects these defects may be the actual cause 

of the lattice curvature in this system. Other interest- 
ing systems in which almost perfect conformal crystals 
have been generated include ferrofluid foams in magnetic 
fields [ir] and soap foams ^8] ; in some cases the resulting 
structure contains internal defects. It is hoped that the 
present work may give some insight into their role. 
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